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Abstract 

The  Picayune  Strand  Restoration  Project  is  one  of  many  components  of  the 
Comprehensive  Everglades  Restoration  Project  (CERP)  intended  to  restore 
nearly  700  hectares  of  a  failed  residential  development  in  southwestern 
Collier  County,  FL,  to  its  predevelopment  wetland  conditions.  A  detailed 
analysis  was  performed  to  derive  a  restoration  plan  that  will  achieve  this 
goal.  As  required  by  the  Water  Resources  Development  Act  (WRDA)  2000, 
the  U.S.  Army  Corps  of  Engineers  (USACE)  is  required  to  ensure  that  no 
component  of  CERP  results  in  an  effective  taking  of  land  by  adversely 
impacting  the  level  of  flood  protection  of  adjacent  landowners.  To  ensure 
the  current  level  of  flood  protection  is  maintained,  a  hydrologic  model  was 
developed  to  assess  the  potential  for  flooding  and  to  refine  the  proposed 
flood  mitigation  features.  The  USACE  physically  based  Gridded  Surface 
Subsurface  Hydrologic  Analysis  (GSSHA)  model  was  selected  for  this  effort. 
The  GSSHA  model  simulates  fully  coupled  rainfall  distribution,  extraction, 
retention,  overland  flow,  and  one-dimensional  channel  flow.  Models  of 
varying  resolution  were  developed  from  existing  and  proposed  design  data 
and  were  initially  populated  with  parameter  values  from  a  previous 
hydrodynamic  modeling  effort.  Parameters  were  then  tuned  to  observed 
stage  and  flow  data  using  the  Secant  Levenberg-Marquardt  method,  a 
nonlinear  least  squares  minimization  computer-based  local  search  method. 
The  calibrated  model  is  capable  of  reproducing  canal  flows,  canal  stages, 
and  overland  stages  with  very  high  Nash  Sutcliffe  Forecast  Efficiencies, 
generally  0.9  or  higher.  Subsequent  uncertainty  analysis  allowed  water 
stages  to  be  estimated  with  95%  certainty.  Modeling  and  uncertainty 
analysis  results  allowed  for  refinement  of  the  proposed  flood  mitigation 
features.  The  hydrologic  models  and  analysis  demonstrated  that  some  of  the 
features  in  the  original  plan  were  either  unnecessary  or  overdesigned  and 
could  be  modified  or  eliminated,  resulting  in  $40M  in  flood  control  feature 
construction  cost  savings. 
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1  Introduction 

Purpose 

The  purpose  of  this  report  is  to  describe  the  model  development,  simula¬ 
tions,  and  results  of  Gridded  Surface  Subsurface  Hydrologic  Analysis 
(GSSHA)  numerical  modeling  in  support  of  level-of-service  flood  analysis 
for  refinement  of  flood  protection  features  proposed  in  the  Project 
Implementation  Report  (PIR)  for  the  Picayune  Strand  Restoration  Project 
(PSRP). 

Background 

The  PSRP  will  functionally  restore  more  than  700  hectares  of  land  in 
southwestern  Collier  County,  FL,  to  its  predevelopment  condition  as  part  of 
the  Comprehensive  Everglades  Restoration  Plan.  Formerly  known  as  the 
Southern  Golden  Gate  Estate,  the  project  area  was  intended  as  a  residential 
subdivision  in  the  1950s.  Roads  and  drainage  canals  were  constructed  in  the 
1960s  and  early  1970s  for  that  purpose.  The  planned  residential  develop¬ 
ment  did  not  materialize,  and  the  roads  and  four  canals  overdrain  the  area. 
This  overdrainage  reduces  aquifer  recharge,  increases  freshwater  discharges 
to  downstream  estuaries,  promotes  invasive  upland  vegetation,  interrupts 
ecological  connectivity,  degrades  habitat,  and  increases  the  frequency  of 
wild  fires.  The  recommended  plan  for  the  PSRP  is  to  neutralize  the  existing 
infrastructure  of  the  subdivision  and  restore  its  predrainage  hydrology  and 
ecological  function  for  beneficial  effects  on  flora  and  fauna  in  the  project 
area  and  surrounding  public  lands.  The  structural  features  of  the  PSRP  are 
the  construction  of  a  series  of  pump  stations,  levees,  spreader  berms,  and 
canal  plugs  to  impede  water  flowing  through  existing  canals  and 
redistribute  it  across  the  landscape.  The  U.S.  Army  Corps  of  Engineers 
(USACE),  Jacksonville  District  (SAJ),  in  partnership  with  South  Florida 
Water  Management  District  (SFWMD),  is  constructing  these  features. 
Engineering  support  is  required  for  hydrologic  and  hydraulic  design  of  the 
project  features.  The  PSRP-recommended  plan  was  selected  by  the  Project 
Delivery  Team  as  having  the  greatest  ecological  benefit  at  the  highest  cost- 
benefit  ratio  (USACE  2004).  The  recommended  plan  was  selected  from 
numerous  alternative  plans  after  years  of  rigorous  study. 
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With  the  plan  selected,  the  USACE  is  required  to  maintain  the  existing 
level  of  flood  protection  to  potentially  affected  persons/areas  once  the 
project  is  built.  In  addition,  all  Federal  projects  must  also  be  designed  with 
a  100-year  level  of  flood  protection.  This  report  describes  the  analysis  of 
the  flood  protection  features  only.  Previous  reports  describe  the  ecological 
analysis  and  process  of  selecting  the  preferred  alternative  (USACE  2004). 

In  order  to  adequately  determine  the  design  parameters  for  the  flood 
mitigation  features,  a  tiered- model,  design-analysis  approach  was  applied. 
A  large-scale,  lower-resolution  (Tier  1)  model  enabled  analysis  of  the 
system  at  a  regional  level  and  provided  the  ability  to  estimate  the  response 
of  the  system  in  terms  of  flood  stages  and  flows  relative  to  historical  and 
synthetic  storm  events.  This  Tier  1  (Ti)  model  was  also  used  to  provide 
boundary  conditions  for  more  detailed,  high-resolution  inset  (Tier  2) 
models.  The  Tier  2  (T2)  models  enabled  detailed  analysis  of  flow  and  stage 
behavior  in  the  vicinity  of  the  flood  mitigation  features  proposed  in  the 
recommended  plan.  The  purpose  of  the  T2  models  was  to  determine 
whether  adverse  water  stages  are  anticipated  in  and  around  private  lands 
for  various  return  intervals  (i.e.,  the  5-,  10-,  25-,  and  100-year  events)  and 
to  determine  if  the  proposed  mitigation  features  meet  the  flood  protection 
requirements.  If  protection  features,  such  as  levees,  were  found  to  be 
necessary,  then  the  design  criteria  were  further  refined  at  this  level, 
including  the  evaluation  of  risk  and  uncertainty. 

The  GSSHA  spatially  distributed  hydrologic  model  (Downer  et  al.  2005) 
developed  and  supported  at  the  U.S.  Army  Engineer  Research  and 
Development  Center,  Coastal  and  Hydraulics  Laboratory  (ERDC-CHL), 
was  selected  to  represent  the  PSRP  existing  and  proposed  project 
conditions  at  the  Ti  and  T2  scales.  GSSHA  provides  comprehensive 
treatment  of  watershed  hydrology,  the  capacity  to  substantively  reflect  the 
previously  mentioned  structural  features,  and  an  ability  to  extract  state 
information  from  a  coarse  scale  model  to  serve  as  boundary  conditions  for 
a  finer  scale  inset  model. 

Hydrologic  simulator 

GSSHA  is  a  physics-based,  spatially  explicit  hydrologic  model  with  the 
capacity  to  simulate,  on  a  continuous  basis,  numerous  processes  relevant  to 
the  hydrologic  response  of  a  watershed  system  subject  to  meteorological 
forcing,  including,  among  possible  others,  rainfall  distribution,  plant 
interception,  surface  retention,  evapotranspiration,  vertical  infiltration, 
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two-dimensional  (2D)  overland  flow,  one-dimensional  (lD)  channel  flow, 

2D  groundwater  flow  for  the  unconfined  aquifer  system  and  related  surface 
water-groundwater  interaction,  lake/reservoir  levels,  and  snow  accumula¬ 
tion  and  melt.  Processes  are  treated  with  varying  degrees  of  fidelity  with 
related  computational  and  input  data  requirements,  allowing  the  modeler  to 
choose  a  process  solution  method  that  is  consistent  with  project  require¬ 
ments  and/or  to  explore  the  potential  tradeoffs  between  simulation 
accuracy  and  related  resource  requirements.  Spatial  data  products  such  as 
digital  elevation  models,  surveyed  channel  cross  sections,  soil  classification, 
vegetative  cover,  land  use,  and  hydro-geological  characterizations  of  the 
subsurface  can  readily  be  incorporated  to  support  model  development  and 
model  process  parameterization. 

The  GSSHA  model  structure  is  also  designed  to  allow  the  modeler  to 
reflect  into  the  model,  in  a  site-specific  manner,  various  man-made 
alterations  to  the  hydrologic  and  hydraulic  watershed  system  including, 
among  others,  groundwater  pumping,  subsurface  storm  and  tile  networks, 
culverts,  weirs,  roads,  and  levees.  Other  types  of  external  boundary 
conditions,  such  as  inflow  hydrographs,  or  the  treatment  of  inset  models 
whose  boundary  conditions  are  obtained  from  a  coarser-scale  GSSHA 
model,  are  also  possible.  The  model  also  has  the  capability  to  actively 
simulate  the  exchange  between  the  channel  system  and  the  floodplain. 
GSSHA-simulated  hydrologic  states  can  be  written  to  files  during  program 
execution  for  subsequent  display  and  evaluation.  (The  interested  reader  is 
directed  to  the  GSSHA  wiki  (http://www.gssahwiki.com)  for  a  complete  summary 
of  the  GSSHA  simulation  model  capabilities,  including  tutorial  documents 
with  related  data  sets,  and  links  to  obtain  a  current  copy  of  the  model 
executable,  which  is  in  the  public  domain  and  supported  on  multiple 
computing  platforms.) 

Interface 

The  Watershed  Modeling  System  (WMS)  is  a  graphically  based,  spatial 
data  processing  environment  that  supports  the  complete  GSSHA  model 
deployment  process.  It  greatly  enhances  the  usability  of  the  GSSHA  model 
for  the  practicing  hydrologic  and  hydraulic  engineer.  Project-specific 
spatial  data  products  that  characterize  the  spatial  distribution  of  physical 
characteristics  of  the  watershed  or  relevant  system  attributes  can  be 
readily  imported  into  the  WMS  and  subsequently  transformed  into 
GSSHA  model  specific  features.  These  features  can  be  displayed  and  easily 
parameterized  or  incorporated  into  simulated  model  processes.  Geo- 
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referenced  imagery  can  also  be  imported  into  WMS  and  overlaid  with 
supplementary  project-specific  spatial  data  sets,  possibly  representative  of 
predetermined  design/impact  features.  Coupling  project  information  with 
the  graphical  display  of  relevant  simulated  states  permits  the  high-level 
evaluation  of  project  alternatives.  The  GSSHA  model  structure  with  its 
interface  within  WMS  enables  a  more  transparent  simulation  and 
evaluation  of  project  alternatives  than  is  possible  with  traditional 
hydrologic  simulation  models,  viz.,  lumped  and  semidistributed  model 
structures.  The  WMS  was  used  extensively  not  only  to  develop  but  also 
evaluate  the  PSRP  Ti  and  T2  GSSHA  models  representing  existing  and 
proposed  project  conditions. 
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2  Preferred  Restoration  Alternative 

Currently  the  project  area  is  drained  with  a  network  of  canals.  A  network  of 
surface  streets  intended  to  allow  access  to  residences  disrupts  overland 
flow.  In  the  Tentatively  Select  Plan,  the  canals  are  plugged,  and  most  of  the 
roads  are  degraded  (Figure  l).  The  flow  conveyance  currently  provided  by 
the  canals  is  replaced  with  pumping  stations  and  spreader  berms,  intended 
to  restore  sheet  flow  (Figure  1).  Restoration  features  include  the  following: 

1.  The  Prairie  Canal  is  plugged  in  several  locations  along  the  length  of  the 
canal. 

2.  Merritt,  Faka  Union,  and  Miller  canals  remain  open  upstream  of  pump 
stations  constructed  as  part  of  the  plan  and  are  plugged  in  several 
locations  downstream  of  the  pump  stations. 

3.  A  small  section  of  the  Faka  Union  Canal  is  left  open  near  the  Port  of  the 
Isles  (POI). 

4.  Pump  stations  are  constructed  on  Merritt,  Faka  Union,  and  Miller  canals. 
These  pump  stations  discharge  water  from  the  upstream  canal  to  the 
spreader  berms,  which  discharge  water  by  gravity  through  a  series  of  weirs 
in  the  berm  creating  sheet-flow  downstream  of  the  berms  for  wetland 
restoration  purposes. 

5.  Most  surface  streets  are  degraded. 

As  shown  in  Figure  1,  the  project  includes  three  areas  requiring  flood 
control:  the  Private  Lands,  the  6Ls,  and  the  POI.  Design  flood  protection 
features  include  the  following: 

1.  The  tie-back  levee  alignment  starts  east  of  the  Merritt  Canal  and  stops 
1,829  m  west  of  the  Miller  Pump  Station. 

2.  A  new  interior  drainage  canal  is  to  provide  drainage  for  the  Private  Lands 
area  north  of  the  tieback  levee  by  conveying  storm  water  from  a  point 
where  the  existing  drainage  canal  starts  at  62nd  Ave  SE  to  the  Miller  pump 
station  intake  canal,  where  it  will  be  pumped  south  into  the  PSRP. 

3.  A  levee  is  built  to  protect  the  6Ls  area. 
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Figure  1.  Preferred  alternative. 
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3  Modeling  Approach 

The  modeling  approach  was  to  (l)  develop  the  Existing  Condition 
(Condition  l)  Ti  model,  (2)  calibrate  and  verify  that  model  to  observed 
stage  and  discharge  data  from  two  historical  tropical  storm  events  that  hit 
the  area,  (3)  include  the  features  for  the  preferred  alternative  analysis 
(Conditions  2),  (4)  build  the  T2  models  (Conditions  1  and  2)  that  derive 
boundary  conditions  from  the  larger  Ti  model,  and  (5)  use  this  suite  of 
models  to  simulate  a  series  of  return  period  storms  events. 

The  Ti  and  T2  Existing  Conditions  (Condition  1)  GSSHA  model  domains 
are  shown  in  Figure  2.  The  figure  also  shows  other  features  including 
relevant  canals  and  structures,  U.S.  Interstate  75  (I-75),  U.S.  Highway 
(Hwy)  41,  and  the  location  of  culverts  along  the  U.S.  Hwy  41  road  way.  At 
the  Ti  level,  the  Condition  1  (Cl)  model  is  the  preproject  condition  without 
any  of  the  PSRP  features  in  place.  All  of  the  canals  are  intact;  no  pumps 
are  installed;  no  flood  protection  features  are  in  place. 

The  Proposed  Project  Model  Condition  model  (Condition  2)  is  shown 
schematically  in  Figure  3.  The  Condition  2  (C2)  model  includes  the 
preferred  alternative  features  described  in  Section  2.  As  shown  in  the 
Figure,  many  of  the  canals  in  Figure  2  are  plugged  in  Figure  3.  These 
canals  are  not  included  in  the  model  as  lD  stream  features  in  the  C2 
model.  Water  from  the  remaining  canals  is  conveyed  via  a  network  of 
pumping  stations  and  spreader  berms  (Figure  3).  In  Figure  3,  the  pumping 
stations  are  located  on  the  canals  behind  the  spreader  berms. 

To  simplify  the  following  discussion,  the  models  will  be  referred  to  first  as 
tier  level  (Ti  or  T2),  followed  by  condition  (Cl)  or  (C2)  (e.g.,  T1C2). 

Tier  1  (TI)  model 

The  Ti  model  provides  the  basis  for  the  more  refined  T2  models.  It 
includes  sufficient  detail  to  provide  large  scale  impacts  and  boundary 
conditions  for  the  T2  models.  Important  features  of  the  Ti  model  are 
described  below. 
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Figure  2.  T1  and  T2  model  domains  and  existing  (Cl)  features. 
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Figure  3.  T  1  and  T  2  model  domains  and  proposed  (C2)  model  features. 


CSUEarms. 


BR37 


■BRQrtjgfJ 

Itneiisianas] 


ALLIGATOR  ALLEY/l-75 


Rrivate 


[lands 


Legend 


— |  Tier  1  Model 
—I  Domain 

— l  Private  Lands  Tier  2 
- -  Model  Domain 

POI  Tier  2  Model 
■»!  Domain 

- 1  6LTier  2  Model 

- ,!  Domain 

Existing  6L  Farms 
Levee 


i  Separated 

Canals 

Major  Roads 

Calibration/ 
Verification  Stations 

Other  US-41 
Bridges 


Project  Features 

Spreader  Berms 

Tie  Back  Levee 

Degraded  SGGE 
Roads 

■  Proposed  6L  Levee 


■  Canal  Plugs 
[0  Pump  Station 


Domain  and  boundary  conditions 

The  domain  for  the  Ti  model  was  selected  based  on  available  data,  the 
need  for  coherent  boundary  conditions,  and  inclusion  of  important  areas. 
The  northern  boundary  of  the  model  coincides  with  I-75  (Figure  2).  The 
interstate  has  an  embankment  and  a  drainage  canal  with  structures  that 
control  flow  into  the  model  domain.  A  no-overland  flow  condition  was 
prescribed  at  the  northern  boundary.  Flow  into  the  domain  from  the  area 
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upstream  of  the  model  domain  is  only  through  the  canals  in  the  existing 
conditions  and  into  the  remnant  canals  immediately  upstream  the  pump 
stations  in  C2.  These  boundary  flows  were  extracted  from  a  previous 
hydrodynamic  model  applied  for  the  selection  of  the  preferred  alternative 
(USACE  2004).  Flows  into  the  canals  were  input  as  inflow  hydrographs 
into  the  lD  stream  network  for  Cl  and  as  overland  flow  in  C2  models. 

The  southern  overland  flow  boundary  of  the  model  is  a  specified  head 
boundary  representative  of  the  tidal  stage  fluctuations  over  the  period  of 
simulation,  and  it  was  derived  from  the  model  applied  to  select  the 
preferred  alternative  (USACE  2004).  The  eastern  and  western  boundaries 
are  no-overland  flow  conditions,  supported  by  local  topography.  Topo¬ 
graphy  data  were  derived  from  a  number  of  sources  and  composited  into  a 
series  of  topographic  maps  at  various  grid  scales  of  3, 5, 10, 50, 100  m,  etc. 
WMS  was  used  to  determine  the  location  of  the  topographic  breaks.  Model 
grid  cell  elevations  were  based  on  these  same  topographic  data.  Total  area 
in  the  computational  domain  is  928  km2. 

Grid  size  and  time-step 

A  grid  size  analysis  was  performed  to  select  the  optimal  grid  size  for  the  Ti 
model.  In  GSSHA,  the  optimal  grid  size  captures  all  the  relevant  spatial 
data,  without  undue  computational  burden.  A  limitation  of  the  grid  size 
was  the  need  to  process  and  simulate  the  data  on  desktop  computers.  A 
number  of  grid  sizes,  from  30  to  200  m,  were  tested.  The  final  grid  size  of 
the  Ti  model  was  120  m.  The  simulation  of  the  model  with  this  grid  size 
was  facilitated  by  the  inclusion  of  the  overland  flow  structures  that 
represent  the  bridges/culverts  under  U.S.  Hwy  41,  which  are  on  the  order 
of  30  m  in  length.  The  total  number  of  grid  elements  in  the  domain  is 
64,920.  With  the  overland  flow  structures  included  in  the  model,  this  grid 
size  was  more  than  adequate  to  capture  all  salient  spatial  variability  in  the 
model  domain.  With  this  grid  size,  the  model  runs  with  time-steps  on  the 
order  of  10  s.  The  development  of  the  overland  flow  structures  in  GSSHA 
is  discussed  in  Section  4. 

Hydraulic  structures  and  controls 

The  PSRP  Cl  models  contain  numerous  hydraulic  structures  and  other 
overland  flow  control  features,  as  described  above.  In  GSSHA,  the  varying 
structures  are  simulated  in  the  following  ways: 
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•  Canals:  All  canals  were  simulated  explicitly  in  the  model  as  lD  channel 
flow.  Cross-section  and  bottom  elevation  data  were  available  from  the 
as-built  drawings  and  defined  the  canal  cross  sections  and  bottom 
elevations  in  the  models.  The  same  time  varying  stage  applied  at  the 
southern  overland  flow  boundary  is  used  to  specify  the  stage  at  the 
model  outlet,  the  Faka  Union  Canal.  This  canal  extends  to  the  model 
boundary.  Canals  received  specified  flow  hydrographs  from  the 
northern  boundary  and  lateral  flow  from  the  overland.  Flow  exchanges 
between  the  canal  and  the  overland  cells  occur  in  proportion  to  the 
water  stage  in  the  canal  and  surrounding  overland  cells.  When  the 
stage  in  the  canal  exceeds  the  ground  surface  elevation  and  the  top  of 
bank  of  the  canal,  water  can  flow  out  of  bank  onto  the  overland  flow 
plane,  depending  on  the  water  surface  elevation  on  the  overland  flow 
plane.  When  the  stage  in  surrounding  overland  cells  is  higher  than  the 
bank  elevation  of  the  channel  and  if  the  stage  in  the  channel  is  lower 
than  that  of  the  surrounding  cells,  overland  water  flows  into  the  canal. 
Weirs  in  the  canals  were  included  in  the  model.  Relevant  physical 
information  was  obtained  from  the  previous  hydrodynamic  model 
(USACE  2004).  Discharge  coefficients  were  initially  prescribed  from 
this  model  and  then  calibrated  to  observed  data.  Canals  that  were 
either  plugged  or  physically  removed  as  part  of  the  preferred 
alternative  were  removed  from  the  stream  network. 

•  Roads.  Levees.  Spreader  Canals,  and  other  Embankments:  Overland 
flow  impediments,  described  above,  are  included  in  the  GSSHA  model 
as  embankments.  Embankments  block  water  across  faces  on  the 
overland  flow  grid.  Embankments  block  overland  flow  if  the  water  level 
on  the  overland  flow  plane  does  not  exceed  the  embankment  height  in 
the  face  of  the  cell  blocking  the  water  flow.  When  the  embankment 
height  is  exceeded,  flow  over  the  embankment  is  computed  by  a  broad- 
crested  weir  equation.  Water  surface  elevations  on  both  sides  of  the 
embankment  are  analyzed  for  backwater  effects  during  these 
computations.  Hydraulic  structures  through  embankments  are 
represented  as  overland  flow  structures  with  specified  rating  curves. 
The  embankment  for  U.S.  Hwy  41  contained  37  of  these  overland  flow 
structures.  Elevations  for  existing  structures  were  obtained  from  as- 
built  drawings.  Elevations  for  proposed  structures  were  obtained  from 
design  drawings.  As  a  result  of  this  modeling,  some  of  the  structure 
sizes,  locations,  and  elevations  were  modified  and  reflected  back  into 
the  modeling  efforts. 
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•  Overland  Roughness  and  Retention:  Overland  roughness  is  specified 
based  on  the  1995  land  use  and  vegetation  map.  The  model  contains 
18  land  use/vegetation  types.  These  land  uses  are  employed  to  assign 
overland  hydraulic  roughness  coefficients  (Manning’s  coefficient)  for  the 
model.  The  primary  land  use/vegetation  types  in  the  model  are  hydric 
flatland  (26%),  swamp  forest  (15%),  cypress  (12%),  mangrove  (12%),  wet 
prairie  (9%),  and  unknown  (10%).  Roughness  coefficients  were 
originally  assigned  based  on  the  previous  hydrodynamic  model  applied 
during  the  alternative  selection  process  and  then  calibrated  to  historic 
events,  as  described  below.  As  described  by  the  predominant  vegetation, 
the  PSRP  domain  is  largely  a  swamp,  and  the  area  has  depression 
features  at  both  the  grid  and  subgrid  scale.  These  depression  features 
retain  water  at  low  overland  flow  stages  (io-2  m  magnitude)  but  have 
little  effect  on  flood  flows,  which  may  be  several  meters  in  depth.  During 
model  development,  it  was  decided  to  minimize  the  explicit  modeling  of 
the  fine-scale  depression  features  by  smoothing  out  the  overland  flow 
grid  on  a  macroscale  and  then  including  depression  storage,  specified 
with  an  index  map,  for  lost  depression  storage  during  the  grid 
smoothing  process.  A  detailed  analysis  was  conducted  to  calculate  the 
difference  in  storage  between  the  smoothed  120  m  model  grid  and  the 
original  3  m  topographic  digital  elevation  model  used  to  develop  the 
GSSHA  grid.  The  calculated  depression  storage  was  then  applied  to  the 
GSSHA  model  as  a  map  of  retention  depth.  There  were  42  separate 
depression  storage  areas  identified  and  added  to  the  GSSHA  simulation. 
In  GSSHA,  flow  on  the  overland  flow  plane  will  not  occur  until  the 
retention  depth  has  been  exceeded;  then  flow  above  the  retention  depth 
is  calculated  with  the  Manning  formula. 

•  Soils  and  Infiltration:  While  the  event-based  simulations  of  the  suite  of 
return  period  storms  are  simulated  with  a  high  water  table  assumption 
(i.e.,  no  infiltration),  the  need  to  establish  initial  conditions  for 
calibration  and  verification  required  that  soils  and  infiltration  be 
considered  in  the  model.  Soils  were  derived  from  the  1990  Natural 
Resources  Conservation  Service  survey  for  Collier  County,  FL,  as 
provided  by  the  SFWMD.  Four  different  soil  classifications  were  used  in 
the  model:  Immokalee,  Boca,  Riviera,  and  Sanibel. 

Tier  2  (T2)  models 

T2  models  were  developed  to  refine  and  represent  with  higher  fidelity 
some  of  the  relevant  hydraulic  features  in  the  proposed  alternative.  They 
were  applied  to  refine  the  hydraulic  design  criteria  for  the  proposed  flood 
protection  features  and  to  analyze  potential  adverse  effects  of  the  areas 
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where  protection  is  needed.  Three  T2  inset  models  were  developed  for  the 
existing  (Cl)  and  proposed  project  conditions  (C2):  Private  lands,  6Ls,  and 
the  POI.  The  model  domains  and  relevant  features  are  shown  in  Figures  4 
through  8,  with  Figures  4,  6,  and  8  showing  the  Cl  features  for  the  Private 
Lands,  6Ls,  and  POI  models,  respectively,  and  Figures  5,  7,  and  8  showing 
the  C2  features. 


Figure  4.  T1C1  Model  domain  for  the  Private  Lands. 
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Figure  5.  T2C2  model  domain  for  the  Private  Lands. 
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Figure  6.  T2C1  domain  and  features  for  the  6Ls. 
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Figure  7.  T2C2  domain  and  features  for  the  6Ls. 
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Figure  8.  T2C1  and  T2C2  domain  and  features  for  the  POI. 
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The  T2  models  share  many  common  properties.  The  domains  and  grid 
sizes  for  the  models  were  selected  to  allow  the  features  in  the  respective 
models  to  be  simulated  accurately  and  efficiently.  Specified  head  overland 
flow  boundary  conditions  derived  from  the  Ti  models  were  specified  along 
all  the  boundaries  for  all  the  models,  unless  otherwise  noted  below.  All 
three  models  were  constructed  with  a  20  m  grid  size  and  time-steps  on  the 
order  of  10  s.  Specifics  for  each  model  are  listed  in  Table  1.  All  simulation 
times  refer  to  a  single  event  simulation  on  an  Intel  Xeon  X5760  CPU  at 
2.93  GHz. 
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Table  1.  Tier  2  models  properties. 


T2  Model 

Model  Area 
(km1 2 3) 

T1  Area 

(%) 

No.  of  Grid 
Cells 

Simulation  Time 
(hr) 

Private  Lands 

67 

7 

167,648 

4 

6Ls 

155 

17 

388,731 

6.5 

Port  of  the  Isles 

41 

4 

102,492 

2.5 

Many  data  sets  in  the  T2  models  were  imported  from  the  Ti  model.  Among 
these  were  soil  properties,  land  use,  roughness  coefficients,  retention 
depths,  and  any  other  parameter  whose  representation  in  these  two 
domains  was  not  dependent  on  model  resolution.  In  the  Cl  models,  the  grid 
resolution  allowed  features  such  as  roads  to  be  explicitly  included  in  the 
overland  flow  model.  In  C2  models,  these  features  were  assumed  to  be 
degraded.  If  a  canal  traversed  the  model  domain,  that  canal  was  included  in 
the  T2  model.  In  all  such  cases,  the  upstream  boundary  condition  was  an 
input  hydrograph,  and  the  downstream  boundary  condition  was  a  time- 
varying  specified  head.  Items  specific  to  each  of  the  T2  models  are  described 
below. 

Private  Lands  model 

The  T2  model  domain  for  the  Private  Lands  area  is  shown  in  Figures  4  and 
5  for  Cl  and  C2,  respectively.  The  areas  that  cannot  be  adversely  impacted 
by  the  proposed  project  are  highlighted  in  the  figures.  Features  specific  to 
this  model  include  the  following: 

1.  The  overland  flow  boundary  condition  on  the  southern  boundary  in  C2 
model  was  modified  to  a  constant  flow  boundary  condition  to  remove 
water  from  the  domain  instead  of  creating  what  was  considered  an 
artificial  backwater  effect  due  to  overland  stages  imposed  at  this  boundary. 

2.  In  C2,  the  Miller  Canal  pump  system  was  implicitly  simulated  by  a  time 
series  of  outflows  assigned  to  model  cells.  This  flow  time  series  was 
representative  of  how  the  pump  system  would  operate  during  the  storm 
events  simulated  in  this  work  (USACE  2004).  The  outflow  time  series  for 
the  Miller  pump  system  are  shown  in  Figure  9. 

3.  In  C2,  the  proposed  alternative  includes  the  tieback  embankment  and  the 
spreader  berm  as  depicted  in  Figure  5.  The  spreader  berm  includes  a  series 
of  10  weirs  designed  to  convey  water  from  inside  the  internal  spreader 
berm  area  into  the  landscape  in  a  sheet  flow  regime.  T2  model  resolution 
allowed  for  the  inclusion  of  these  hydraulic  details  in  the  model.  The 
design  document  drawings  for  the  spreader  berm,  embankments,  and 
weirs  were  used  to  define  these  features. 
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Figure  9.  Discharge  time  series  for  the  Miller  pump  system. 


6Ls  model 

The  model  domain  for  the  T2  6Ls  area  is  shown  in  Figures  6  and  7. 

Enclosed  within  the  6Ls  levee  are  farm  lands  that  shall  not  be  adversely 

impacted  by  the  PSRP  restoration  features  and  its  operations.  Features 

specific  to  this  model  include  the  following: 

1.  The  only  hydraulic  structures  applicable  to  Cl  and  C2  are  the  U.S.  Hwy  41 
culverts  (Figures  6  and  7).  There  is  one  additional  culvert  included  in  the 
C2  model  in  the  southern-most  segment  of  the  USACE  proposed 
protective  levee  (Figure  7). 

2.  Embankments  in  Cl  include  the  U.S.  Hwy  41  road  and  the  6Ls  levee. 

3.  In  C2,  the  proposed  alternative  includes  a  protective  levee  (Figure  7)  that  is 
offset  approximately  60  m  from  the  6Ls  levee  alignment.  The  goal  of  this 
feature  is  to  provide  flood  protection  to  the  areas  enclosed  by  the  6Ls  levee. 

Port  of  the  Islands  (POI)  model 

The  domain  for  the  POI  T2C1  and  T2C2  models  is  shown  in  Figure  8. 

Features  specific  to  this  model  include  the  following: 
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1.  For  C2,  there  is  a  no  flow  boundary  condition  at  the  upstream  end  of  the 
canal  since  in  the  preferred  alternative,  all  canals  segments  south  of  the 
pump  stations  and  north  of  this  junction  will  be  plugged.  For  the  most 
downstream  node  of  the  canal,  a  stage  boundary  condition  was  assigned 
consistent  with  the  Ti  model. 

2.  The  hydraulic  structures  in  Cl  and  C2  are  the  U.S.  Hwy  41  culverts 
(Figure  8)  and  an  existing  weir  in  the  Faka  Union  Canal  located  just  north 
of  U.S.  Hwy  41  (not  shown).  The  proposed  alternative  does  not  consider 
additional  hydraulic  or  control  features  in  this  area  beyond  the  existing 
conditions. 

3.  The  portion  of  the  Faka  Union  Canal  that  will  remain  open  in  the  preferred 
alternative  is  the  only  channel  within  the  model  domain.  This  applies  to 
both  conditions  since  the  T2  model  domain  bisects  the  Faka  Union  Canal 
at  the  junction  of  the  three  main  canals.  The  difference  in  how  the  canal  is 
treated  in  Cl  and  C2  resides  in  the  boundary  conditions  as  described 
previously.  No  other  canal  is  explicitly  simulated  in  these  conditions. 

4.  The  preferred  alternative  does  not  alter  the  existing  topography  in  the 
vicinity  of  POI;  therefore,  topography  in  Cl  and  C2  was  represented  with 
only  one  source  of  data. 

5.  The  only  embankment  in  Cl  and  C2  is  U.S.  Hwy  41.  In  the  PIR  (USACE 
2004),  a  protective  embankment  feature  was  being  considered  in  the  area 
around  POI.  T2C2  model  simulations  showed  that  this  proposed  feature 
did  not  significantly  add  flood  protection  beyond  the  existing  conditions. 
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4  Gridded  Surface  Subsurface  Hydrologic 
Analysis  (GSSHA)  Model  Enhancements 
for  Picayune  Strand  Restoration  Project 
(PSRP) 


The  PSRP  presented  many  unique  challenges  to  the  modeling  team.  The 
project  required  models  at  varying  levels  of  size  and  resolution,  where  the 
finer-resolution  T2  models  are  subsets  of  the  larger,  coarser  Ti  model 
domain.  To  function  properly,  the  higher-resolution  T2  models  had  to 
extract  temporally  and  spatially  varying  overland  and  channel  boundary 
conditions  from  the  larger  Ti  models.  To  facilitate  the  development  of 
boundary  conditions  from  the  Ti  model  for  application  in  the  T2  models, 
the  GSSHA  model  was  modified  such  that  when  the  Ti  model  was  executed, 
the  Ti  model  output  the  required  spatially  and  temporally  varying 
information  needed  for  the  T2  model  simulations.  The  T2  model  was  then 
executed  using  the  output  from  the  Ti  model  as  an  input  for  boundary 
conditions.  The  Ti  model  writes  output,  overland  depth  and  canal  flow,  for 
Ti  model  cells  located  in  proximity  to  the  T2  model  boundary.  The  location 
of  the  T2  model  boundary  is  a  required  input  for  the  Ti  model  to  produce 
the  correct  output.  The  T2  model  reads  the  time  varying  overland  depths 
and  canal  flows  from  the  Ti  model  output  file  and  uses  them  for  the 
overland  head  and  canal  flow  boundaries.  Because  the  grid  resolution  of  the 
two  models  is  different,  the  overland  depth  values  from  the  Ti  model  must 
be  interpolated  to  the  cells  on  the  boundary  of  the  T2  model  grid.  An 
inverse-distance-squared  interpolation  method  is  employed. 

Another  important  development  was  the  inclusion  of  point  flow  structures 
on  the  overland  flow  plane.  A  major  highway,  U.S.  Hwy  41,  traverses  the 
lower  portion  of  the  model  domain  (Figures  1  and  2).  Because  the  highway 
traverses  a  flood  prone  area,  it  is  elevated  approximately  3  m.  Several 
bridges/box  culverts  exist  along  the  road  embankment  (Figures  1  and  2). 
Compared  to  the  resolution  of  the  Ti  models,  these  features  are  relatively 
small  and  are  not  well  represented  as  breaks  in  the  embankment  that 
represents  U.S.  Hwy  41  in  the  Ti  models.  Refining  the  Ti  model  to  capture 
these  features  in  the  overland  flow  grid  was  not  computationally  reasonable. 
To  allow  the  proper  simulation  of  these  important  features,  the  ability  to 
simulate  point  overland  flow  control  structures  was  added  to  the  GSSHA 
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model.  The  overland  point  structures  function  by  passing  water  from  one 
overland  flow  cell  to  another  based  on  a  rating  curve  that  relates  either 
depth  or  difference  in  water  surface  elevation  between  the  cells  to  flow. 
Overland  cells  that  contain  structures  are  specified  on  the  overland  flow 
boundary  condition  map  and  identified  with  a  unique  number  that  is 
associated  with  a  rating  curve.  These  overland  point  structures  represented 
flow  through  the  hydraulic  structures  along  the  U.S.  Hwy  41  road 
embankment  (Figure  1).  The  rating  curves  applied  were  developed  from  the 
MIKE  SHE/11  hydrologic/hydraulic  model  used  for  the  alternative 
assessment  modeling  (USACE  2004). 
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5  PSRP  Model  Calibration 

Thirteen  parameters  were  specified  as  adjustable,  and  their  names, 
meanings,  and  prescribed  lower  and  upper  bounds  are  listed  in  Table  2. 
The  first  parameter  listed  in  Table  2,  viz.,  ch_n,  holds  the  uniform 
Manning’s  roughness  n  value  for  all  channel  links  reflected  in  the  TiCi 
model.  The  second  parameter  listed  in  Table  2,  viz.,  dcf,  holds  the  uniform 
value  for  the  forward  discharge  coefficient  for  all  five  weirs  that  are 
reflected  in  the  TiCi  model.  As  the  canals  have  been  in  place  for  50  years 
without  maintenance,  the  stream  bed  elevation  values  derived  from  the  as- 
built  drawings  have  considerable  uncertainty  and  were  adjusted  with  the 
third  parameter  listed  in  Table  2,  viz.,  adj,  which  is  the  value  to  uniformly 
adjust  all  of  the  link  node  elevation  values  specified  in  the  channel  input 
file  for  the  T1C1  model.  Prior  manual  model  calibration  efforts  with  the 
T1C1  model  resulted  in  the  specification  of  this  adjustable  model 
parameter.  The  fourth  parameter  listed  in  Table  2,  viz.,  adj2,  holds  the 
value  to  uniformly  adjust  bankfull  depth  values,  from  their  specified  base 
values,  for  each  channel  link  represented  in  the  T1C1  model.  This 
parameter  was  specified  to  be  adjustable  based  on  expert  solicitation1.  The 
embankments  along  the  canals  are  the  result  of  spoil  from  the  canals  and 
are  not  uniform.  The  bankfull  depth  represents  an  effective  bankfull  depth. 
Base  overland  roughness  values  were  uniformly  modified  via  the  value  for 
the  specified  adjustable  model  parameter  named  multi.  Base  retention 
depth  values  were  uniformly  modified  via  the  value  for  the  specified 
adjustable  model  parameter  named  mult2.  Based  on  experience  obtained 
via  prior  manual  model  calibration  efforts  with  the  T1C1  model,  the  value 
for  the  specified  adjustable  model  parameter  named  rain_mult  permitted 
all  input  precipitation  values  to  be  increased,  to  a  small  reasonable  degree, 
in  attempts  to  improve  model  to  measurement  misfit  during  calibration. 
The  values  associated  with  the  remaining  final  six  specified  adjustable 
model  parameters,  viz.,  rtc_mult_36,  rtc_mult_37,  rtc_mult_52, 
rtc_mult_55,  rtc_mult_w,  and  rtc_mult_e  pertain  to  uniform  adjustment 
of  their  respective  base  rating  curve  values  for  the  stage-discharge 
relationships  for  the  culverts  along  U.S.  Hwy  41,  as  indicated  in  Table  2. 


1  Charles  W.  Downer,  2012,  personal  communication,  ERDC/CHL,  Vicksburg,  MS. 
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Table  2.  Specified  adjustable  parameters  for  the  T1C1  model,  their  description,  and  specified 

lower  and  upper  bounds. 


Parameter 

Parameter  description 

Lower 

Bound 

Upper 

Bound 

ch_n 

Channel  Manning's  n  value 

0.025 

0.075 

dcf 

Forward  discharge  coefficient  for  weirs  in  channel 
network 

1 

5 

adj 

Adjustment  to  all  node  elevation  values  in  the 
channel  input  file 

-0.2 

0.2 

adj2 

Adjustment  to  bankfull  depth  values  in  the  channel 
input  file 

-1 

1 

multi 

Multiplier  for  original  overland  roughness  values 

0.2 

5 

mult2 

Multiplier  for  original  retention  depth  values 

1.00E-01 

1.00E+01 

rain_mult 

Multiplier  for  original  input  rainfall  values 

1 

1.1 

rtc_mult_36 

Multiplier  for  original  rating  curve  values  for  36 

1.00E-02 

1.00E+02 

rtc_mult_37 

Multiplier  for  original  rating  curve  values  for  37 

1.00E-02 

1.00E+02 

rtc_mult_52 

Multiplier  for  original  rating  curve  values  for  52 

1.00E-02 

1.00E+02 

rtc_mult_55 

Multiplier  for  original  rating  curve  values  for  55 

1.00E-02 

1.00E+02 

rtc_mult_w 

Multiplier  for  remaining  original  rating  curve  values 
west  of  FU* 

1.00E-02 

1.00E+02 

rtc_mult_e 

Multiplier  for  remaining  original  rating  curve  values 
east  of  FU 

1.00E-02 

1.00E+02 

*Faka  Union 


Estimation  of  the  13  specified  adjustable  parameters  for  the  T1C1  model 
was  undertaken  by  matching  observed  and  simulated  flows  and  stages  at 
Faka  Union  1  (FU-i),  and  also  by  matching  observed  and  simulated  stages 
at  four  ditch  gaging  locations  along  U.S.  Hwy  41,  viz.,  US4i-Bridge36, 
US4i-Bridge37,  US4i-Bridge52,  and  US4i-Bridge55,  for  a  single  storm 
event  (Tropical  Storm  Jerry:  20  Aug  1995, 14:00  -  26  Aug  1995, 11:00), 
resulting  in  830  observations  for  the  calibration  process.  Locations  of  the 
calibration/ verification  stations  FU-i,  US4i-Bridge36,  US4i-Bridge37, 
US4i-Bridge52,  and  US4i-Bridge55  are  shown  in  Figure  2.  Six  groups 
consisting  of  observations  and  their  model-simulated  counterparts 
composed  the  objective  function  as  indicated  in  Table  3.  The  objective 
function,  which  is  the  quantitative  measure  of  model  to  measurement 
misfit,  is  defined  in  the  weighted-least-squares  sense.  Weight  assignment 
for  each  observation  group  was  adjusted  such  that  at  the  beginning  of  the 
computer-based  calibration,  each  observation  group  was  seen  to  be  of 
equal  importance. 
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Table  3.  Unique  observation  groups  defined  to  calibrate  T1C1  model. 


Observation  Group  Name 

Description 

Number  of  Values 

mf 

Modeled  flows  at  FU-1 

142 

ms 

Modeled  stages  at  FU-1 

142 

ms36_l 

Modeled  stages  at  US41-Bridge36 

135 

ms37_l 

Modeled  stages  at  US41-Bridge37 

134 

ms52_l 

Modeled  stages  at  US41-Bridge52 

138 

ms55_l 

Modeled  stages  at  US41-Bridge55 

139 

The  Secant  Levenberg-Marquardt  (SLM)  method,  a  nonlinear  least 
squares  minimization  computer-based  local  search  method  (Skahill  et  al. 
2009),  as  implemented  by  Skahill  et  al.  (2012),  was  employed  to  efficiently 
calibrate  the  T1C1  model.  With  the  exception  of  the  two  parameters  named 
adj  and  adj2,  which  can  take  on  negative  values,  the  remaining  13 
specified  adjustable  model  parameters  were  log  transformed  in  attempts 
to  better  accommodate  scaling  issues  resulting  from  the  use  of  different 
units  for  different  parameters  and  to  decrease  the  degree  of  nonlinearity  of 
the  parameter  estimation  problem;  past  experience  has  demonstrated  that 
greater  efficiency  and  stability  of  the  parameter  estimation  process  can 
often  be  achieved  through  this  means  (Skahill  et  al.  2009). 

The  SLM  local  search  with  the  T1C1  model  was  configured  for  efficiency.  In 
particular,  no  full  update  of  the  model  sensitivity  matrix  was  employed 
together  with  single  column  cyclic  updating1.  The  SLM  local  search 
completed  after  33  forward  model  calls,  reducing  the  objective  function 
from  an  original  value  of  34,582  to  a  final  value  of  22,141.  The  final 
estimated  model  parameter  set  is  listed  in  Table  4  together  with  an  estimate 
of  its  uncertainty.  Note  that  the  confidence  limits  provide  only  an  indication 
of  parameter  uncertainty.  They  rely  on  a  linearity  assumption  which  may 
not  extend  as  far  in  parameter  space  as  the  confidence  limits  themselves. 
Clearly,  the  rating  curve  multipliers,  viz.,  the  specified  adjustable  model 
parameters  named  rtc_mult_36,  rtc_mult_3y,  rtc_mult_52,  rtc  mult  55, 
rtc_mult_iv,  and  rtc_mult_e  were  identified  with  the  most  uncertainty 
relative  to  the  remaining  model  parameters.  Computed  model  parameter 
sensitivities  to  all  observations  are  listed  in  Table  5  (Doherty  2004).  Model- 
simulated  results  are  most  sensitive  to  the  parameters  named  adj  and  adj2. 
The  information  shown  in  Figure  10  and  listed  in  Table  6  summarizes 


1  Refer  to  Skahill  et  al.  (2009)  and  Skahill  et  al.  (2012)  for  details  regarding  the  SLM  method  and  its 
implementation. 
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model-to-measurement  misfit  at  the  final  estimated  parameter  set.  The 
computed  Nash  and  Sutcliffe  efficiency  scores  indicate  a  T1C1  calibrated 
model  of  predictive  value. 


Table  4.  Final  estimated  parameter  set  for  the  T1C1  model  with  uncertainty  estimate. 


Parameter 

Estimated  Value 

95%  Confidence  Limits 

Lower  Limit 

Upper  Limit 

ch_n 

0.045844 

0.044531 

0.047195 

dcf 

1.61301 

1.521033 

1.71055 

adj 

0.001531 

-0.004223 

0.007285 

adj2 

0.000042 

0.00002 

0.000063 

multi 

1.100165 

0.988629 

1.224284 

mult2 

0.856535 

0.83239 

0.881381 

rain_mult 

1.065394 

1.043216 

1.088043 

rtc_mult_36 

0.93399 

0.682793 

1.277602 

rtc_mult_37 

0.952339 

0.608473 

1.490534 

rtc_mult_52 

0.901433 

0.698453 

1.163402 

rtc_mult_55 

0.965863 

0.763623 

1.221665 

rtc_mult_w 

0.999786 

0.684533 

1.460226 

rtc_mult_e 

0.96641 

0.699074 

1.335978 

Table  5.  Computed  model 
parameter  sensitivities  to  all 
observations. 


Parameter 

Sensitivity 

ch_n 

1.409143 

dcf 

2.01872 

adj 

47.02873 

adj2 

624.9205 

multi 

0.741138 

mult2 

1.830177 

rain_mult 

3.064292 

rtc_mult_36 

0.62694 

rtc_mult_37 

0.87159 

rtc_mult_52 

1.406658 

rtc_mult_55 

1.438024 

rtc_mult_w 

1.302043 

rtc_mult_e 

0.537431 
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Figure  10.  Comparison  of  all  simulated  and  observed  data  (FU-1  flow  and  stage,  stage 
values  at  ditch  locations  numbered  36, 37,  52,  and  55)  for  the  calibrated  T1C1  model. 


Table  6.  Summary  statistics  of  model  to  measurement  misfit  at  the  final  estimated  parameter  set 
for  the  T1C1  model  across  different  observation  groups,  including  their  subdivision  or  aggregation. 


Data  Comparisons 

R2 

Nash-Sutcliffe 

Efficiency 

Scores 

Mean 

Error 

(m) 

Mean 

Absolute 

Error 

(m) 

RMSE2 

(m) 

FU-1  discharge  data1 

0.98 

0.97 

1.121 

3.931 

4.831 

FU-l  stage  data 

0.97 

0.97 

0.00 

0.03 

0.04 

All  ditch  stage  data 

0.90 

0.87 

-0.02 

0.07 

0.11 

Ditch  numbered  36  stages 

0.86 

0.77 

-0.07 

0.10 

0.14 

Ditch  numbered  37  stages 

0.86 

0.77 

-0.02 

0.12 

0.15 

Ditch  numbered  52  stages 

0.98 

0.97 

0.02 

0.04 

0.05 

Ditch  numbered  55  stages 

0.99 

0.99 

0.00 

0.02 

0.03 

discharge  errors  in  m3  s1. 
2Root-mean-square  error. 
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6  PSRP  Model  Predictions 

Hydrologic  models  are  typically  calibrated  by  adjusting  parameters 
encapsulated  in  the  simulator  until  there  is  an  acceptable  level  of 
agreement  between  a  set  of  historical  data  and  their  model  simulated 
counterparts.  The  parameters  obtained  via  calibration  are  often  used  in 
the  model  to  predict  system  behavior  for  one  or  more  predefined  scenarios 
of  interest.  Regardless  of  the  calibration  method  employed  and  the  type  of 
model  applied  (e.g.,  empirical  or  physics  based),  some  if  not  all  of  the 
parameter  values  obtained  through  the  calibration  process  possess  a 
degree  of  quantifiable  uncertainty.  This  is  because  the  observed  data 
contain  measurement  errors  and  also  because  the  model  never  perfectly 
represents  the  watershed  system  or  exactly  fits  the  observation  data. 
Where  model  parameters  are  uncertain,  so  too  will  be  model  predictions. 

For  the  PSRP  alternatives  modeled  at  the  Ti  level,  nonlinear  calibration 
constrained  predictive  analysis  was  performed  (Cooley  and  Vecchia  1987; 
Vecchia  and  Cooley  1987).  Due  to  computational  constraints,  project 
alternatives  evaluated  at  the  T2  level  only  permitted  a  less-intensive  linear 
uncertainty  analysis  (Doherty  2004).  However,  due  to  time  limitations, 
linear  uncertainty  analyses  were  also  performed  with  the  T1C2  model  for 
six  locations  along  the  proposed  tieback  levee  system.  These  uncertainty 
analyses  that  were  performed  for  the  PSRP  project  alternatives  supported 
risk-based  hydrologic  design  in  that  the  upper  95%  prediction  limit  was 
computed  for  stages  at  specific  locations  of  interest  in  each  case.  Levee 
superiority  analysis  was  performed  to  assess  whether  the  levee  crest  was  at 
least  0.6  m  greater  than  the  95%  confidence  level  for  the  100-year  event. 

Nonlinear  calibration  constrained  predictive  uncertainty  analysis  (Cooley 
and  Vecchia  1987;  Vecchia  and  Cooley  1987)  involves  maximizing/ 
minimizing  a  specific  prediction  of  interest,  s,  while  keeping  the  calibration 
model  effectively  in  a  calibrated  state.  In  particular,  the  objective  function 
value  associated  with  the  estimated  parameter  set,  o(p) ,  is  relaxed  to 

accommodate  for  the  desired  level  of  confidence,  1  -  a ,  as  indicated  in  the 
following: 


0(p)-  <£>(p)<  m<y1rFa{m,n  -m) 


(1) 
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®(p).  m,  cr2 ,  Fa(m,n-m),  and  n  represent  the  objective  function  for  a 

parameter  set  p,  the  number  of  specified  adjustable  model  parameters,  the 
reference  variance,  the  F  distribution,  and  the  number  of  non-zero 
weighted  observations  employed  to  calibrate  the  model.  The  simultaneous 
prediction  confidence  interval  is  given  by 


Prob 


s  -  y[mFa {m,  n  -  m)as  <s<s  +  y[mFa (m,n  -  m)cr. 


=1  —  a 


(2) 


where  <js  is  the  standard  deviation  of  s.1  Application  of  this  methodology 

generally  involves  the  development  of  a  composite  model  that  is 
comprised  of  the  model  run  under  calibration  conditions  followed  by  the 
predictive  model  representative  of  a  given  alternative  scenario  run  under 
predictive  conditions.  Simulation  of  the  predictive  model  yields  the 
specific  unique  prediction  that  is  to  be  maximized/minimized,  as 
mentioned  above. 

While  tractable  for  many  model  deployments,  nonlinear  predictive  analysis 
is  nonetheless  nontrivial  in  terms  of  its  model  simulation  time  require¬ 
ments.  If  forward  model  execution  times  are  prohibitively  expensive,  then 
one  can  employ  a  linear  analysis  wherein  the  only  additional  postcalibration 
computational  burden  is  to  compute  the  sensitivity  of  the  prediction  of 
interest  to  the  set  of  specified  adjustable  model  parameters.  This  will 
require  at  a  minimum  (m  +  l)  forward  model  calls.  The  variance  of  the 
uncertainty  of  the  prediction  in  this  case  is  given  by 

o j  =  yt{XtC~1{e)X  +  C-1(p))  1y  (3) 

where  y,  X,  C  (e),  and  C  (p)  represent  the  sensitivity  of  the  prediction  to  the 
specified  adjustable  model  parameters,  the  model  sensitivity  matrix  at  the 
estimated  parameter  set,  the  covariance  of  the  error  term  for  measurement 
noise,  and  the  covariance  matrix  associated  with  the  estimated  parameter 
set.  All  of  the  input  in  the  equation  immediately  above,  with  the  exception  of 
the  prediction  sensitivity  vector  y,  is  information  that  was  obtained  at  the 
end  of  the  computer-based  model  calibration  for  the  TiCi  model. 

Nonlinear  calibration  constrained  predictive  uncertainty  analysis 
simulations  were  performed  to  compute  the  maximum  simulated  stage,  at 


1  For  complete  details,  see  Cooley  and  Vecchia  (1987)  and  Vecchia  and  Cooley  (1987). 
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the  95%  confidence  level,  at  two  grid  cell  locations  for  the  T1C2  POI  model. 
The  noted,  two  GSSHA  model  grid  cell  locations  are  indicated  in  Figure  11. 
The  predictive  model  was  forced  with  precipitation  data  representative  of 
the  5-day,  100-year  design  condition.  Table  7  summarizes  the  two  nonlinear 
predictive  analysis  runs. 

Figure  11.  The  T1C1  uncertainty  analysis  assessment  points. 


Table  7.  Maximum  predicted  stage  values,  at  the  95% 
confidence  level,  simulated  by  the  T1C2  GSSHA  model 
forced  with  the  5-day,  100-year  design  storm. 


Predictive  GSSHA  Model 
Grid  Cell  Location  (i,j) 

Maximum  Predicted  Stage 
(m) 

171,  179 

1.62 

184,  179 

1.64 

Five  linear  predictive  uncertainty  analyses  were  performed  to  compute  the 
maximum  predicted  stage,  at  the  95%  confidence  level,  at  five  distinct 
locations  associated  with  the  T2  level  GSSHA  Private  Lands  model.  These 
locations  are  designated  in  Figure  12.  As  with  the  POI  predictive  model, 
the  composite  predictive  Private  Lands  model  was  forced  with  precipita¬ 
tion  data  representative  of  the  5-day,  100-year  design  storm.  For  each  of 
the  five  locations,  the  vector  y  in  Equation  (3)  above  had  to  be  determined 
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from  the  composite  predictive  Private  Lands  model.  The  computation  of 
the  prediction  sensitivity  vector  y  in  each  case  required  six  forward  model 
calls  of  the  composite  predictive  Private  Lands  model.  Only  six  forward 
model  calls  were  required  to  populate  the  vector  y  because  the  specified 
adjustable  model  parameters  named  dcf,  adj2,  rtc_mult_36,  rtc_mult_3y, 
rtc_mult_52,  rtc_mult_55,  rtc_mult_w,  and  rtc_mult_e  were  all  held 
fixed  at  their  estimated  values  as  listed  in  Table  4  because  by  virtue  of  the 
project  alternative  they  were  known  to  not  be  relevant.  At  each  location, 
after  the  vector  y  was  determined  via  first-order  finite  difference  calcula¬ 
tions,  Equation  (3)  was  subsequently  used  to  compute  the  uncertainty  of 
the  prediction  of  maximum  stage  at  the  95%  confidence  level.  Table  8 
summarizes  the  linear  uncertainty  analyses  for  the  five  locations  for  the 
composite  predictive  Private  Lands  model. 


Figure  12.  T2C2  Private  Lands  model.  The  black  dots  indicate  the  approximate  locations 
where  linear  predictive  uncertainty  analyses  were  performed  to  compute  the  maximum 
simulated  stage  at  the  95%  confidence  level  when  the  model  is  forced  with  the  design  storm. 
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Table  8.  Maximum  predicted  stage  values,  at  the  95%  confidence 
level,  simulated  by  the  T2C2  Private  Lands  GSSHA  model  forced 
with  the  5-day,  100-year  design  storm. 


Predictive  T2  GSSHA  Private  Lands 
Model  Grid  Cell  Location  (i,j) 

Maximum  Predicted  Stage 
(m) 

153,  311 

3.40 

153,  370 

3.46 

136,  179 

3.18 

152,  200 

3.14 

152,  269 

3.29 

Six  linear  predictive  uncertainty  analyses  were  performed  to  compute  the 
maximum  predicted  stage,  at  the  95%  confidence  level,  at  six  distinct 
locations  associated  with  the  T2  level  GSSHA  6Ls  model.  These  locations 
are  designated  in  Figure  13.  As  with  the  previous  predictive  models,  the 
6Ls  model  was  forced  with  precipitation  data  representative  of  the  5-day, 
100-year  design  storm.  For  each  of  the  six  locations,  the  vector  y  in 
Equation  (3)  above  had  to  be  determined  from  the  composite  predictive 
6Ls  model.  The  computation  of  the  prediction  sensitivity  vector  y  in  each 
case  required  14  forward  model  calls  of  the  composite  predictive  6Ls 
model.  At  each  location,  after  the  vector  y  was  determined  via  first-order 
finite  difference  calculations,  Equation  (3)  was  subsequently  used  to 
compute  the  uncertainty  of  the  prediction  of  maximum  stage  at  the  95% 
confidence  level.  Table  9  summarizes  the  linear  uncertainty  analyses  for 
the  six  locations  for  the  composite  predictive  6Ls  model. 

A  representation  of  the  T1C2  GSSHA  model  is  depicted  in  Figure  14.  Six 
linear  predictive  uncertainty  analyses  along  with  a  single  nonlinear  analysis 
were  performed  to  compute  the  maximum  predicted  stage,  at  the  95% 
confidence  level,  at  six  distinct  locations  associated  with  the  T1C2  model. 
These  locations  are  shown  in  Figure  14,  and  the  maximum  water  elevation 
in  the  corresponding  grid  cells  is  listed  in  Table  10.  The  predictive  T1C2 
model  was  forced  with  precipitation  data  representative  of  the  5-day,  100- 
year  design  storm.  For  each  of  the  six  locations,  for  the  linear  analysis,  the 
vector  y  in  Equation  3  above  had  to  be  determined  from  the  predictive 
model.  The  computation  of  the  prediction  sensitivity  vector  y  in  each  case 
required  14  forward  model  calls  of  the  predictive  model.  At  each  location, 
after  the  vector  y  was  determined  via  first-order  finite  difference  calcula¬ 
tions,  Equation  (3)  was  subsequently  used  to  compute  the  uncertainty  of  the 
prediction  of  maximum  stage  at  the  95%  confidence  level.  Table  10 
summarizes  the  linear  and  nonlinear  uncertainty  analyses  for  the  six 
locations  for  the  T1C2  GSSHA  model. 
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Figure  13.  T2C2  6Ls  model.  The  black  dots  indicate  the  approximate  locations  where  linear 
predictive  uncertainty  analyses  were  performed  to  compute  the  maximum  simulated  stage  at 
the  95%  confidence  level  when  the  model  is  forced  with  the  design  storm.  Two  separate 
analyses  were  performed  at  distinct  grid  cell  locations  in  the  vicinity  of  the  southernmost 

black  dot  indicated  in  the  figure. 
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Table  9.  Maximum  predicted  stage  values,  at  the  95% 
confidence  level,  simulated  by  the  T2  predictive  6Ls  model 
forced  with  the  5-day  100-year  design  storm. 


Predictive  T2  GSSHA  Miller 
Model  Grid  Cell  Location  (i,j) 

Maximum  Predicted  Stage 
(m) 

379,  279 

2.438 

354,  215 

2.028 

461,  260 

2.183 

582,  174 

4.002 

411,  279 

2.300 

579, 175 

2.470 
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Figure  14.  T1C2  model.  The  six  pink  dots  indicate  the  approximate  locations  where  linear  and 
nonlinear  predictive  uncertainty  analyses  were  performed  to  compute  the  maximum  simulated 
stage  at  the  95%  confidence  level  when  the  T1C2  model  is  forced  with  the  design  storm. 


Table  10.  Maximum  predicted  stage  values,  at  the  95% 
confidence  level,  simulated  by  the  T1  level  GSSHA  predictive 
C2  model  forced  with  the  5-day  100-year  design  storm. 


Predictive  T2  GSSHA  Miller 
Model  Grid  Cell  Location  (i,j) 

Maximum  Predicted  Stage 
(m) 

24,  205 

3.685  (linear),  3.692 
(nonlinear) 

33,  188 

3.635 

41,  126 

3.199 

43,  148 

3.279 

43,  181 

3.320 

111,  114 

2.392 
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7  Maximum  Flood  Levels  and  Implications 
to  the  Protection  Features  Design 

To  quantify  the  residual  flooding  caused  by  the  proposed  PSRP  restoration 
features,  the  T2  values  of  gridded  maximum  flood  stage  in  Cl  were 
subtracted  from  the  T2C2  flood  stages  for  all  storm  event  frequencies 
analyzed  in  this  study.  The  results  of  this  comparison  along  with  the  PSRP 
model  predictive  analysis  discussed  in  the  previous  section  allowed  for 
further  refinement  of  the  proposed  restoration  features. 

The  cell-by-cell  difference  of  the  maximum  flood  extent  between  T2C2  and 
T2C1  for  the  100-year  storm  event  is  shown  at  POI  in  Figure  15.  As  shown 
in  the  figure,  this  analysis  indicated  either  no  change  in  water  level  or  a 
lower  water  level  for  the  C2  model  relative  to  the  Cl  model.  These  results 
indicate  that  the  proposed  protection  feature  that  was  conceptualized  as  a 
protective  levee  just  east  of  the  developed  areas  north  of  U.S.  Hwy  41  does 
not  provide  any  additional  flood  protection  relative  to  the  existing 
conditions.  Moreover,  the  proposed  restoration  features  for  the  PSRP  will 
most  likely  alleviate  flooding  conditions  in  this  vicinity.  Results  for  all 
other  storm  events  showed  similar  trends  to  the  100-year  frequency  storm. 
A  water  surface  profile  for  the  10-year  and  100-year  storm  events  just  east 
of  the  developed  areas  showed  that  water  levels  for  C2  are  consistently 
lower  than  Cl  (Figure  16). 

Doing  a  similar  analysis  for  the  Private  Lands  T2C2  model  showed  that  the 
tie-back  levee  protective  feature  that  was  originally  conceptualized  in  the 
PIR  (Figure  17)  extending  approximately  3,000  m  west  of  the  Miller  pump 
station  to  the  south  of  the  Private  Lands  area  could  be  reduced  to 
approximately  1,800  m  in  length  without  decreasing  the  level  of  flood 
protection  this  feature  was  designed  to  provide.  This  is  illustrated  in 
Figure  17,  where  the  difference  in  maximum  flood  extent  between  C2  and 
Ci  is  shown  on  a  cell-by-cell  basis.  These  results  show  that  increases  in 
water  depths  in  C2  relative  to  Cl  do  not  extend  to  the  private  lands  areas. 
Water  surface  profiles  along  the  recommended  levee  extension,  Figure  17, 
are  shown  in  Figure  18. 
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Figure  15.  Cell-by-cell  difference  of  maximum  flood 
extent  for  the  POI  T2  model. 


Figure  16.  Water  surface  profiles  along  the  east  boundary  of  the  developed  areas  in  the  POI 

T2  model. 
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Figure  17.  Cell-by-cell  difference  of  maximum  flood  extent  for  the  Private  Lands  T2 
model.  Western  extent  of  the  protective  tie-back  levee  system  as  proposed  in  the 
PIR  and  as  recommended  in  this  study  also  shown. 
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A  water  surface  profile  just  south,  along  the  alignment  of  the  western 
portion  of  the  tie-back  levee  system,  shows  that  one  of  the  reasons  why  the 
water  levels  do  not  increase  in  C2  relative  to  Cl  near  the  private  lands 
areas  is  due  to  the  existence  of  a  spoil  mound  located  approximately 
540  m  east  from  the  western  end  of  the  system.  This  spoil  mound  runs 
approximately  1000  m  south  and  prevents  water  from  the  spreader  berm 
from  flanking  the  proposed  levee  where  it  could  reach  the  Private  Lands 
area.  The  resolution  selected  for  these  T2  models  captured  all  these 
important  details  allowing  for  a  refinement  of  the  protective  feature  design 
as  originally  conceptualized  in  the  PIR. 

The  cell-by-cell  maximum  flood  extent  difference  grid  for  the  6Ls  T2 
model  is  shown  in  Figure  19.  Due  to  the  way  in  which  water  is  spread  over 
the  landscape  by  the  proposed  restoration  features  (i.e.,  the  Miller  pump 
station  and  spreader  berm  system),  the  flooding  extent  in  areas  just 
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outside  the  6Ls  farm  are  shown  to  increase  in  C2  relative  to  Cl.  The 
residual  flooding  east  of  the  proposed  6Ls  protective  feature  varies  from 
0.15  m  east  of  the  farms  to  1.0  m  in  areas  just  south  of  the  restoration 
features.  A  water  surface  profile  (Figure  20)  along  the  proposed  feature  for 
the  10-year  and  100-year  storm  event  frequencies  also  validates  the  results 
in  Figure  19.  The  proposed  embankment  is  therefore  necessary  to  protect 
the  6Ls  farms  from  flooding.  Additional  refinement  of  flood  control 
features  for  this  area  may  be  needed  to  alleviate  backwater  effects  or 
seepage  into  the  6Ls. 

Based  on  the  maximum  water  levels,  the  Private  Lands  levee  was  modified 
as  noted,  and  the  POI  levee  is  not  being  constructed.  These  project 
modifications  result  in  approximately  $40  million  (U.S.  dollars)  savings  in 
project  construction  costs  without  a  reduction  in  flood  protection. 

Figure  18.  Water  surface  profiles  along  the  western  portion  of  the  proposed  tie-back  levee  system 

(see  Figure  17  for  profile  locations). 
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Figure  19.  Cell-by-cell  difference  of  the  maximum  flood  extent  between  C2  and  Cl  in 

the  6Ls  T2  model. 
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Figure  20.  Water  surface  profiles  along  the  proposed  6Ls  protective  levee.  Profile  starts  at  the 
northern  end  of  the  feature  and  ends  at  the  intersection  of  the  feature  with  U.S.  Hwy  41. 
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8  Summary 

The  GSSHA  hydrologic  model  was  applied  to  execute  a  multitiered,  level- 
of-service  flood  study  analysis  on  the  existing  PSRP  and  multiple 
idealizations  of  the  preferred  alternative  for  the  flood  protection  features. 
The  GSSHA  models  were  built  from  existing  spatial  data  and  design 
documents.  The  models  included  a  number  of  hydraulic  structures, 
including  canals,  weirs,  roads,  levees,  spreader  canals,  and  pumps.  The 
simulation  of  flows  through  embankments  on  the  overland  flow 
computational  domain  was  facilitated  by  the  development  of  overland 
hydraulic  structures  feature  in  GSSHA  that  allows  flow  to  be  passed  from 
one  overland  flow  grid  cell  to  the  next  based  on  rating  curves  of  depth  or 
difference  in  stage  versus  flow.  The  model  also  utilized  a  new  capability  in 
the  GSSHA  model  to  derive  boundary  conditions  for  inset  models  and 
boundary  conditions  from  the  coarser  Ti  model  to  drive  the  finer- 
resolution  T2  models.  Initial  model  parameters  were  derived  from 
previous  large-scale  coarser  hydrodynamic  models  used  in  the  analysis  to 
select  the  preferred  alternative.  The  models  were  calibrated  to  existing 
flow  and  stage  data  in  the  canals  and  overland  flow  domains.  Parameters 
specified  as  adjustable  as  part  of  the  T1C1  computer-based  model  calibra¬ 
tion  included  not  only  model  process  related  parameters  but  also 
parameters  related  to  model  configuration,  input  precipitation  forcing, 
and  adjustment  of  the  base  rating  curve  values  that  controlled  the  flow 
through  the  numerous  bridge  crossings  along  U.S.  Hwy  41.  Hence,  the 
postcalibration  model  predictive  uncertainty  analysis  included  these 
elements  in  the  quantification  of  the  uncertainty  associated  with  specific 
stage  predictions  selected  to  be  analyzed  with  each  project  alternative. 

One  of  two  tractable  methodologies  for  postcalibration  model  predictive 
uncertainty  analysis  was  used  with  each  of  the  PSRP  GSSHA  models 
representing  project  alternatives.  Employment  of  either  of  these  two 
analysis  methods  enabled  quantification,  at  the  95%  confidence  level,  of 
the  maximum  simulated  stage  at  specific  predetermined  locations  within 
the  computational  domains  of  the  models  representing  the  proposed 
project  condition.  All  of  the  GSSHA  models  developed  to  characterize 
project  alternatives  were  forced  with  the  5-day,  100-year  design  storm. 
Nonlinear  postcalibration  predictive  uncertainty  analysis  was  employed 
with  the  Ti  project  alternatives.  A  linear  uncertainty  analysis  was 
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employed  to  quantify  the  uncertainty  associated  with  predictions  from  the 
T2  project  alternatives,  primarily  due  to  project  resource  limitations.  The 
results  from  these  analyses  supported  risk-based  hydrologic  design  of 
PSRP  protection  features.  The  hydrologic  models  and  analysis 
demonstrated  that  based  on  maximum  flood  stage  alone,  some  of  these 
features  were  either  unnecessary  or  overdesigned  with  respect  to  the 
original  design  selected  in  the  tentatively  selected  plan.  Suggested 
modifications  could  result  in  approximately  $40  million  savings  in 
construction  costs.  As  this  study  focused  on  maximum  flood  stage  alone, 
any  concerns  about  flooding  duration  and  possible  backwater  or  seepage 
effects  would  need  to  be  addressed  with  additional  analysis. 

This  report  and  the  software  GSSHA  are  products  of  CHL.  For  more 
information  about  GSSHA,  contact  Natalie  Elwart,  Hydrologic  Systems 
Branch,  Coastal  and  Hydraulics  Laboratory,  U.S.  Army  Engineer  Research 
and  Development  Center,  3909  Halls  Ferry  Road,  Vicksburg,  MS  39180- 
6199.  Information  on  GSSHA  can  also  be  found  on  the  GSSHA  wiki  at 
http://gsshawiki.com. 
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